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WHYRAC,  A NEW  MODULAR  ONE-DIMENSIONAL 
EXPLODING  WIRE  CODE 


I.  Introduction 

The  dynamic  evolution  of  an  exploding  wire  or  an  array  of  parallel- 
strung  exploding  wires  depends  in  detail  on  at  least  three  different 
types  of  physical  phenomena:  the  magnet ohydrodynamlc  fluid  motion  of  the 
wire  plasma,  the  coupling  of  the  wire  plasma  to  the  driving  electric 
circuit,  and  the  atomic  physics  and  radiative  energetics  of  the  hlgh-Z 
plasma.  These  different  classes  of  phenomena  are  intimately  cross- 
coupled.  For  example,  the  plasma  temperature/density  profile  depends 
sensitively  on  atomic  excitation  energy  sources  and  sinks,  radiative 
energy  transport,  and  ohmic  energy  deposition  by  the  external  circuit, 
as  well  as  on  MHD  phenomena  such  as  compression  and  shocks.  The  circuit 
current  depends  on  the  plasma  temperature/density  profile,  in  a way  that 
cannot  consistently  be  modeled  in  terms  of  lumped  circuit  elements.  The 
radiative  and  atomic  phenomena,  and  particularly  the  flux,  fluence,  and 
spectrum  of  emitted  radiation,  obviously  depend  strongly  on  the  fluid 
parameters,  i.e.  temperature/density  profile,  of  the  plasma. 

This  report  discusses  a recently  developed  numerical  code,  WHYRAC, 

^^Ich  treats  all  of  these  phenomena  in  a fully  self-consistent  manner. 

WHYRAC  represents  a significant  advance  in  flexibility,  efficiency,  and 

inclusion  of  more  detailed  physics  over  the  first  effort  of  this  type^ . 

WHYRAC  is  one-dimensional,  permitting  variation  of  macroscopic  parameters 

only  in  the  radial  coordinate  in  cylindrical  geometry.  Thus  it  can  treat 

either  a single  wire  (possibly  one  of  the  wires  in  a multiple  wire  array) 
Note:  Manuscript  submitted  February,’  8.  1978. 
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as  a Z-pinch,  i.e.  an  axially  symmetric  plasma  cylinder,  confined  by  an 

azimuthal  magnetic  field  driven  by  an  axial  current,  or  it  can  treat  an 

array  of  wires  mounted  parallel  to  each  other  by  representing  the  array 
as  an  azimuthally  sjmmetric  annulus.  This  report  deals  only  with  the 

physics  and  numerical  techniques  embodied  in  WHYRAC.  Results  from  the 

code  for  specific  physical  problems  and  choices  of  parameters  will  be 

discussed  in  forthcoming  reports. 

WHYRAC  is  a significant  advance  over  all  previous  efforts  of 
this  type.  For  the  first  time,  the  coupling  of  the  wire  plasma  to  the 
external  circuit  is  treated  in  an  exact,  self-consistent  way.  Partially 
as  a result  of  this,  energy  is  conserved  very  accurately  in  the  plasma/ 
generator  system.  The  code  is  modularized  for  extreme  flexibility,  both 
to  handle  a variety  of  situations  and  to  permit  continual  upgrading  of 
the  physics.  For  example,  the  atomic  physics/radiation  package  can  use 
a perfect  gas  or  other  equation  of  state,  coronal  equilibrium,  local 
thermodynamic  equilibrium  (L.T.E.),  or  more  exact  model,  and  presently 
Includes  a frequency-diffusion  radiation  transport  model  which  can  also 
be  even  further  upgraded.  The  external  circuit  can  be  determined  by 
specification  of  either  a driving  voltage  V(t)  or  a circuit  current  I(t). 
The  transport  coefficients  can  easily  be  changed,  in  a way  that  is  con- 
sistent throughout  the  code,  e.g.  various  forms  of  anomalous  resistivity 
can  be  specified,  even  if  they  depend  self-consistently  on  a variety  of 
plasma  and  circuit  characteristics.  Many  of  the  choices  made  in  writing 
the  code,  which  will  be  discussed  below,  can  be  traced  back  to  this 
requirement  of  flexibility. 

The  report  is  divided  into  three  main  parts.  The  first  one  reports 
the  >!HD  equations  to  be  solved  and  some  of  the  problems  associated  with 
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their  solutions.  The  second  part  describes  the  external  circuit  equations 

how  they  are  solved  self -consistently  with  the  MHD  equations.  The  third 
section  discusses  the  atomic  physics  and  radiation  transport  part  of 

the  code.  Although  it  does  not  describe  any  atomic  physics  and  radiation 

package  in  detail,  it  shows  how  the  coupling  is  made  to  the  rest  of  the 

code  and  along  which  lines  the  modularization  takes  place.  A more 

detailed  description  of  the  atomic  physics  code  will  appear  elsewhere. 

II . The  MHD  Equations 

The  MHD  equations  are  well-knoiro  and  can  be  found  in  many  text- 

2 

books.  They  are  presented  here  mainly  for  the  sake  of  completeness. 

The  equations  describe  the  time  evolution  of  a single  fluid,  two-temper- 
ature quasi-neutral  plasma.  The  local  ratio  of  electron  to  ion  density, 
n^/n^,  is  equal  to  the  mean  ionization  state  Z (r,t),  determined  by  the 
atomic  physics  package.  The  continuity  and  momentum  equations  in 
cylindrical  geometry  are 


1°  + i A. 

St  r ;ir 


(rov)  = 0 , 


(1) 


^ — A.  (rpv^)  = - ^ (2) 

^t  r or  Sr  c 


where  q is  the  mass  density,  v is  the  radial  fluid  velocity,  j is  the 

(axial)  current  density,  B is  the  (azimuthal)  magnetic  field,  and  p is 

the  total  pressure  (electrons  + ions). 

WHYRAC  is  a two- temperature  code,  but  the  electron  and  ion  temper- 

tures  T and  T.  are  treated  in  a non-parallel  way  numerically.  The  two 
e i 

primary  quantities  that  are  solved  for  are  T^  and  the  total  plasma  energy 
density  E,  defined  by 
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The  electron  energy  tern  n kT  / (yl)  includes  the  potential  energy  of 
ionization.  The  effective  adiabatic  index  •/,  including  these  contri- 
butions, is  calculated  in  the  atomic  physics  package  and  discussed  in 
Section  IV.  The  two  "energy"  equations  then  read: 


M + i r "(E  + p)  V + q ] = j • g + P, 
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where  q is  the  total  heat  flux,  q = - K 
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^ is  the  longitudinal  electric  field,  = j - Zr^^j-  aCT^-T^),  and 
u is  the  (axial)  current  flow  speed,  defined  by  j = - en^u,  and  P^^^  is 
the  radiative  loss  power.  The  transport  coefficients  and  3^  are 

given  by  Braginskii^,  in  terms  of  various  coefficients  v,  5 ,3  which  he 
tabulates : 

.2 
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Ion  heat  conduction  (generally  small)  was  neglected  in  Eq.  (-J-b) . 

T;  and  p are  treated  as  secondary  quantities.  The  former  is  obtained 
from  Eq.  (3),  while  the  latter  is  obtained  by  using  T^,  T^,  n^  and  n^ 
in  an  equation  of  state  (perfect  gas  is  used),  where  n^  is  calculated 
in  the  atomic  physics  package  described  in  Section  k. 

This  numerical  treatment  affords  great  accuracy  in  overall  energy 
conservation^,  and  accurate  calculation  of  Tg,  which  is  a controlling 
parameter  in  many  sensitive  areas  of  the  code,  e.g.  evolution  of  atomic 
states,  radiation  emission,  resistivity.  T^  is  not  always  calculated 
with  great  accuracy,  but  this  is  not  a significant  deficiency,  since  T^ 
is  only  needed  as  a diagnostic  and  does  not  appear  directly  in  the 
equations  of  evolution. 

The  extra  equation  describing  the  diffusion  of  the  magnetic  field, 
derived  from  the  generalized  Ohm's  law,  is 


M = . (vB)  + . (<a) 

5t  cr  or  r Sr  e Sr  n Sr 

e 


where  T is  the  (cross-field)  resistivity, 

3 

Tl  = 1.15  Z In  A t"^ 

^ ®eV 


whose  classical  value  is 


(6b) 


The  current  density  profile  j(r)  is  then  determined  from  Faraday's  law, 

LA.  rB  (r)  = j rr). 
r Sr  c 
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The  explicit  flux-corrected  transport  (FCT)  algorithm  was  used  to 
solve  Equations  (l)-(5)  numerically,  with  the  exception  of  the  diffusion 
terms.  Those  terms  plus  Eq.  (6)  were  solved  using  an  implicit  tridia- 
gonal solver.  In  the  heat-diffusion  term  q ^ the  thermo-electric  term 
3,  (TbT  ^ ^ removed  and  solved  with  the  convective  terms.  In 

fact,  the  change  in  the  energy  density  due  to  that  term  in  Eq.  (k)  is 
given  by: 

M 4.  i r 3 u = C (?) 

ot  r hr  A 

In  this  equation,  3.  is  of  the  form  3,  = n kT  f (jjd  t ) • Thus,  the 
/.  A e e e e 

thermo-electric  term  looks  like  a convective  term  and  can  be  solved  for 
explicitly  on  the  same  time  scale  as  the  fluid  convective  terms  pro- 
vided that: 

a)  u is  not  larger  than  the  maximum  velocity 
used  in  the  Courant  criterion  (Ivj  + c^) . 

b)  f(j)  T ) is  of  order  1. 

e e 

We  shall  see  that  u does  not  go  much  above  c in  the  course  of  the 

s 

calculation,  where  c is  the  ion  acoustic  speed.  Also,  f (a)  t ) Is 

s e e 

shown  in  Figure  1 and  reaches  a maximum  of  0.55  at  x t =0.12.  The 

e e 

thermo-electric  term  then  does  not  add  any  extra  restrictions  on  the 
fluid  time  step  and  can  be  included  in  an  explicit  calculation  like  the 
one  used  for  the  convective  terms. 

The  equilibration  term  betvmen  electrons  and  ions  in  Eq.  (6)  might 
also  be  a source  of  time-step  problems.  At  the  beginning  of  the  run, 

the  plasma  is  cold  and  the  electron- ion  equilibration  time  can  be  quite 
short,  compared  to  the  time  step  specified  by  other  requirements.  In 
order  not  to  be  limited  by  this  short  equilibration  time,  we  follow  a 
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treatment  similar  to  that  one  used  by  Christiansen  and  Roberts.^  Let 
us  write  this  term  in  the  following  way 


at 


T.-T  m.T  . 

i e 1 ei 

with  T,  = 


(6) 


'1  3(Y-l)m^ 


where  v is  the  adiabatic  index  for  the  electron  thermal  plus  atomic 
excitation  and  ionization  energy,  discussed  in  Sec.  L.  The  equilibration 
time  can  become  tiuch  smaller  than  the  hydrodynamic  time  step  5t  (defined 
by  the  Courant  condition),  and  we  define  AT  = previous 

equation  can  be  written 


at 


with  the  solution 


) 


C9) 


AT  = (AT)^,j 


(1C) 


Going  back  to  the  original  variable  yields 

T = (T  ) (T.)  (1  - (11) 

® ® old  ’■old 

Thus  far,  the  discussion  has  been  quite  general  and  applicable  to 

any  MHD  problem.  At  this  point,  we  discuss  some  of  the  physics  associated 

specifically  with  exploding  wire  and  Z-pinch  problems. 

Since  a region  of  perfect  vacuum  cannot  be  maintained  in  an  Eulerian 

fluid  code,  we  legislate  a minimum  density  for  the  gas  surrounding  the 

Z-pinch,  i.e.  whenever  the  density  falls  below  some  arbitrarily  chosen 

value  in  any  given  cell,  additional  mass  is  introduced  to  raise 

0 to  0 , • Typically  we  choose  o . to  be  1C  of  the  initial  peak 
min  min 

density,  but  other  choices  can  be  made  instead.  This  standard  technique 


does  not  cause  significant  errors  in  mass  conservation,  but  other  types 
of  difficulties  do  occur  in  treating  the  MHD  of  the  low  density  regions 


including  the  following. 

According  to  the  classical  resistivity  formula  (5d),  the  resistivity 

is  essentially  independent  of  density  (for  a fully  ionized  gas),  and  thus 

can  be  very  small  in  a hot,  low  density  region.  If  the  resistivity  is 

small,  the  current  tends  to  flow  at  the  outside  of  the  conducting  region 

(skin  effect).  Thus  a vicious  circle  occurs,  with  the  current  flowing 

exclusively  at  the  outside  of  the  low  density  (o  = d . ) region  (rather 

min 

than  in  the  region  occupied  by  the  Z-pinch  in  a real  experiment),  and 
heating  this  low  density  plasma  to  further  reduce  the  local  resistivity. 

In  real  life,  the  resistivity  in  low  density  regions  becomes 
anomalously  high,  because  of  plasma  instabilities,  thus  allowing  the 

current  to  penetrate  to  higher  density  regions.  Marginal  stability 

8 9 

calculations  ’ have  shown  that  the  exact  value  of  the  anomalous  transport 
coefficients  do  not  control  the  physics  in  this  regime;  rather,  it  is  the 
turn-on  condition  for  instability  to  occur  that  is  crucial.  In  this  spirit, 
we  include  an  anomalous  resistivity  of  the  arbitrarily  chosen  form 


'anom 


m 

'class 


X Max 


'ICC, 


(u/Cs)^"', 


(12) 


2 

whenever  u > c^  = (Z  T j-  T. )/m..  This  prescription  has  the  effect  of 
s e i 1 

limiting  the  electron  drift  speed  u to  close  to  the  ion  sound  speed  c . 

s 

We  OTuld  note,  however,  that  we  believe  that  a much  more  sophisticated 
treatment  of  anomalous  resistivity,  based  on  detailed  study  of  the 
relevant  plasma  instabilities,  will  be  needed  to  explain  some  aspects  of 
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wire  phenomena.  It  is  our  intention  to  pursue  these  aspects  of  the 
physics  and  to  continue  to  upgrade  the  code  here. 


In  addition  to  the  inclusion  of  anomalous  resistivity,  it  is  some- 
times necessary  to  legislate  that  current  may  flow  only  in  regions  where 
p s c^y  where  is  chosen  to  "be  larger  than  by  some  arbitrarily 

chosen  factor.  Whenever  this  is  done,  the  results  are  checked  to  insure 
that  they  do  not  depend  significantly  on  the  choice  of  p^. 

Numerical  difficulties  also  arise  because  the  Alfven  speed  can  become 
very  large  in  the  low  density  region,  as  the  current,  and  therefore  B, 
increases.  As  a result,  the  time  step  (chosen  to  satisfy  the  Courant 

condition)  can  become  intolerably  small.  This  difficulty  is  avoided 

2 2 

by  adding  a term  B /^•p-c^  to  the  inertia  in  the  low-density  region, 

6 7 

thereby  reducing  the  effective  Alfven  speed  ’ . The  quantity  c^^  is 
varied  with  time  to  insure  that  the  Alfven  time  scale  (which  is  really 
an  artifact  of  the  density  floor  in  the  code)  is  no  faster  than  the  true 
hydrodynamic  time  scale.  Of  course,  it  is  important  to  make  sure  that 
this  correction  is  made  only  in  regions  where  the  density  is  lower  than 
the  density  range  of  physical  interest. 


III.  CIRCUIT  EQUATION 

The  external  circuit  (pulse  oower  generator  plus  diode)  is  modeled 

in  the  code  as  a driving  voltage  V (t)  into  an  effective  l umoed  circuit 

8 

with  resistance  Z and  inductance  L„ . The  wire  olasma,  however,  is 
g 8 

treated  in  an  exact  way  as  a distributed  circuit  element,  thus  oermitting 
an  accurate  evaluation  of  the  energy  counled  into  the  plasma.  Our 
circuit  equation  is  then 
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V (t)  = Z I + L 
g g g dt 


P?  ’ 


C13) 


where  V Is  the  voltage  across  the  plasma,  which  is  calculated  by  a 
P 2 

generalization  of  the  method  used  in  Refs.  10.  The  wire  plasma  in  its 
return  current  cage  is  modeled  as  a plasma  cylinder  or  annulus,  in  con- 
tact at  one  end  with  a grounded  perfectly  conducting  metal  return  current 
path,  and  at  the  other  end  with  an  electrode  at  voltage  shown 

in  Fig.  2.  We  then  draw  an  imaginary  current  loop  along  the  plasma  at 
radius  r,  through  the  perfectly  conducting  return  current  path,  and 
across  the  gap  at  Z = as  shown  in  Fig,  2,  and  integrate  Faraday's 
equation 


^ 1 

V X C = - — 

c at 


(1^) 


across  the  loop.  We  then  write  down  the  generalized  Ohm's  law, 

c.  vxB  , , ^ m Z B 

C -i- i 1 (^P  - e V . J-  ? -7  T +6  — X n T - J — 

c n e ® m.  -l  x e |gi  e ^ 


(15) 


where  ^ is  the  electric  field,  p the  electron  pressure,  p.  the  ion 
— e ■ ’ '^1 

pressure  and  the  •?'s  are  defined  as  in  ref,  2. 

From  Eq.  (Ip)  we  note  the  following  two  useful  consequences.  First, 


from  the  axial  component  of  Eq.  (15) 


, J dz  ^jj(r,z)  = ~(r)j(r) 


I - 


, where  ^.  = P,  T is  the  thermo-electric  field, 
c c t e 

Second,  the  radial  component  of  Eq.  (Ip)  shows  that  radial  potential 
drops  ^ 2 dr  ^(r,z)  are  of  the  order  of  within  the  plasma  (i.e.  few 


10 


kilovolts  at  most),  and  therefore  are  very  small  compared  to  axial 
potential  drops  (of  order  few  x IC^  eV  to  few  MeV) , This  is  true  even 
through  the  radial  electric  fields  can  be  large,  if  the  plasma  cylinder 
is  narrow.  It  follows  that  the  current  loop  of  Eq.  (It)  can  be  taken 
at  any  radial  position  within  the  plasma.  The  integrated  form  of 
Eq.  (It)  is  thus 


Vp^  = t(r)j(r)I 


v(r)B(r)£ 


c 


I 


dr'^  B (r;  t). 


(16) 


where  r is  the  wall  radius  and  r is  any  location  within  the  plasma. 
w 

Note  that  the  four  terms  on  the  right  hand  side  of  Eq.  (l6)  play  the 
role  of  resistive,  L,  thermoelectric  and  inductive  circuit  elements,  but 
that  it  is  impossible  to  perform  this  decomposition  into  lumped  circuit 
elements  until  the  current  distribution  is  known,  i.e,  the  circuit 
problem  is  solved. 

It  is  clear  that,  in  general,  the  self-consistent  solution  to 
Eqs.  (6,  15,  16)  can  be  obtained  only  by  iterations  since  I and 

B(r)  are  all  coupled.  V „ is  needed  to  calculate  I from  Eq.  (13)»  the 
radial  profile  of  j or  B is  needed  to  calculate  from  Eq.  (I6),  while 
I is  needed  as  a boundary  condition  in  Eq.  (6)  to  calculate  B (r).  How- 
ever, if  we  take  for  r the  value  corresponding  to  the  cell  location  from 
which  the  current  stops  flowing  by  convention  (and  denoting  this  value 
by  tg) > the  integral  in  Eq.  (I6)  can  be  expressed  in  term  of  I.  If  in 
addition,  we  are  satisfied  with  first-order  accuracy  in  time  (the  overall 
accuracy  of  the  code),  Eq.  (I6)  can  be  written 
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(1^) 


- % 

2 C 


f ” 

J 


where  ^ (r  ) is  the  old  value  of  the  electric  field  at  that  point, 
z c 

Substituting  this  value  in  Eq.  (15),  it  becomes 


fl  + ^ M ^ 


Z I = V 
8 8 


f (r  > 

z 


(16) 


The  iterations  are  not  needed  anymore  because  the  B-field  has  been 
expressed  as  a function  of  I. 

It  remains  to  calculate  the  local  electric  field,  in  order  to 
determine^*  j in  the  plasma  energy  equations  as  mentioned  in  the  pre- 
vious section.  The  most  direct  way  to  evaluate  £ is  given  by  the  axial 
component  of  the  generalized  Ohm's  law,  Eq.  (15),  which  reads. 


<?  = Tl  j-  vB  . -_h 
^ c n 


Sr 


(10) 


However,  use  of  this  form  in  a numerical  calculation  can  lead  to  a 
choppy  electric  field.  For  example,  anomalous  resistivity  can  introduce 
large  differences  from  cell  to  cell.  Also,  low  density  cells  present 
quite  large  variations  in  velocity  and  the  second  term  might,  in  fact, 
lead  to  negative  values  for  the  electric  field. 

Another  expression  for£,  which  is  mathematically  equivalent  but  far 
more  satisfactory  for  numerical  finite  difference  calculations,  can  be 
obtained  directly  from  the  B field.  Eq.  (l6)  can  be  rewritten  as 
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r't 


B dr 


(20 


+ r 


/ 


W 


r 

Since  V is  the  same  every^vhere  in  the  plasma,  ^ inside  the  plasma  can 
PI 

be  calculated  from  point  to  point  starting  with  the  value  given  in 
Eq.  (19) < ^ is  then  obtained  directly  as  a function  of  the  magnetic 

field  flux  and  nothing  else.  B being  the  solution  of  a diffusion 
equation  is  expected  to  be  smooth  and  should  reflect  this  property. 

The  electric  field  obtained  this  way,  when  multiplied  by  j 
(obtained  from  j = t-^  txB),  represents  the  total  electromagnetic  energy 
going  into  the  plasma,  and  should  satisfy  the  conservation  relation 


"J j?  * i " J S 

Eq.  (21)  is  satisfied  to  much  better  accuracy  (of  the  order  of  a few  fS) 

when  £ is  derived  from  the  magnetic  field  flux,  Eq.  (2C),  than  t^en  it 

is  obtained  from  Ohm's  law,  Eq.  (19).  For  the  latter,  the  energy  check 

is  satisfied  only  to  within  5''^  after  ITCC  time  steps.  Also,  in  order 

to  improve  even  further  the  accuracy  of  the  current  determination,  the 

calculation  was  performed  in  two  half  time  steps,  in  order  to  evaluate 

the  contribution  due  to  the  I term 

In  the  introduction,  it  was  mentioned  that  either  V (t)  or  I(t) 

g 

could  be  used  as  input  variables  for  the  driving  circuit.  The  above 
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\ 


>1 


description  applies  to  the  more  complicated  case  when  the  generator 


voltage  is  specified.  For  a current-specified  problem,  the  above  pro- 
cedure simplifies  still  further.  The  plasma  voltage  which  was  calculated 
in  order  to  determine  the  current  does  not  need  to  be  calculated  in  this 
case  (but  still  is  evaluated  for  diagnostic  purposes).  Given  the  current, 
one  can  calculate  the  B field  and  the  current  density.  In  all  cases,  it 
is  found  that 

Jj  atrrdr  = I 

is  satisfied  to  very  good  precision. 


IV.  ATOMIC  PHYSICS  AND  RADIATION 

There  are  two  very  different  motivations  for  computing  the  atomic 
and  radiative  behavior  of  the  wire  plasma  as  a function  of  time.  First 
is  our  obvious  interest  in  the  flux,  fluence,  and  detailed  spectral 
characteristics  of  the  x-radiation  emitted  from  the  plasma.  Second  is 
the  feedback  of  energy  sources  and  sinks  associated  with  atomic  excitation, 
ionization,  and  radiative  transport  on  the  hydrodynamic  (temperature/ 
density)  evolution  of  the  plasma.  Very  different  types  of  atomic/ 
radiative  code  packages  are  required  to  perform  these  two  functions. 

The  detailed  calculation  of  emitted  spectra  is  a diagnostic  function, 
which  may  be  carried  out  quite  infrequently  compared  to  the  time  step  of 
the  code,  (as  a post-processing  facility,  for  example),  but  which 
deserves  very  detailed  (and  therefore  slow  running  and  costly)  treatment 
of  atomic  physics  and  radiative  transport.  In  our  scheme,  this  is 
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accomplished  by  storing  temperature/density  profiles  at  any  chosen  time, 
for  use  as  input  to  separate  atomic/radiation  codes  which  calculate 
detailed  spectra.  On  the  other  hand,  a much  simpler  and  faster  running 
atomic/radiation  package  is  included  in  1</HYRAC  itself  to  accomplish  the 
second  function,  i.e.  to  interface  self-consistently  ^^th  the  energetics 
of  the  MHD  and  circuit  packages  in  the  code.  This  atomic/radiation 
package  must  run  at  every  time  step  of  the  code,  and  thus  depending  on 
its  complexity  might  constitute  a large  fraction  of  the  running  time  of 
the  code.  It  can  be  benchmarked  against  the  more  detailed  diagnostic 
atomic/radiation  codes  to  establish  that  its  accuracy  is  sufficient  unto 
its  needs. 

Our  main  purpose  here  is  to  discuss  the  general  nature  and  inter- 
action of  the  atomic/radiation  section  with  the  rest  of  the  wire  code, 
rather  than  to  discuss  the  specifics  of  these  atomic/radiation  packages. 
This  will  be  done  in  detail  in  subsequent  reports,  as  the  continual  up- 
grading of  accuracy  and  versatility  (e.g.  use  for  different  elements)  of 
these  packages  constitutes  a major  part  of  our  exploding  wire  program. 

We  simply  note  that  the  atomic/radiation  package  presently  being 
used  in  WHYRAC  is  an  equilibrium  model,  in  which  one  excited  state  is 
carried  for  each  ionization  stage,  corresponding  to  the  main  transition 
line  for  that  ionization  level.  Dielectronic  recombination  has  also 
been  included!^  The  treatment  of  radiation  transport  has  been  reduced  to 
a local  one.  For  each  cell  and  for  each  line,  an  opacity  is  calculated 
corresponding  to  the  shortest  optical  path  out  of  the  plasma.  An  escape 
probability  is  then  calculated^^  and  determines  the  fraction  of  the 
radiation  which  escapes  the  system  and  the  fraction  ^^hich  remains  in 
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that  cell.  The  package  needs  as  inputs  only  the  total  ion  density  and 
the  electron  temperature.  The  electron  density,  average  Z,  ionization 
energy  and  radiated  power  (frequency  integrated  or  not)  are  all  calcu- 
lated within  the  atomic/radiation  package. 

Two  problems,  because  of  their  possible  impact  on  the  rest  of  the 
code,  have  received  special  attention.  One  concerns  the  avoidance  of 
negative  temperatures  due  to  radiation  cooling,  the  other  one  the 
avoidance  of  too  low  a value  of  the  adiabatic  coefficient  v at  low 
temperature. 

The  radiation  cooling  term  is  simply 


a 

St 


n kT 
e e 

(v-1) 


rad’ 


(22) 


where  is  evaluated  as  part  of  the  atomic  physics  package.  Eq* 

(22)  can  be  expanded  into 


1 

3/2  n^k 


SE 


chem 


ct 


‘rad 


5/2  n k 
e 


(23) 


5n 

where  r—  is  the  rate  of  change  of  the  electron  density  due  to  the 
dt  ^ 

changes  in  the  ionization  levels,  is  the  potential  energy  of  ioni- 

zation (and,  in  general,  atomic  excitation).  A problem  can  arise  with 
the  last  term,  if  the  change  in  temperature  due  to  that  term  is  to  be 
evaluated  In  the  hydrod3mamic  time  step,  since  the  term  can  become 

large  enough  to  drive  the  temperature  negative,  for  finite  differencing 
on  this  time  scale.  This  would  really  correspond  to  feeding  energy  into 
the  code,  since  the  temperature  is  reset  to  a minimum  positive  value 
whenever  it  goes  negative.  Physically,  of  course,  the  temperature  does 


16 


not  become  negative  because  before  hitting  = 0,  the  radiation  cooling 

term  ^rould  itself  go  to  zero,  so  that  the  temperature  evolution  would 

really  be  like  that  one  shown  in  Figure  Ja.  It  is  not  possible  to  apply 

the  same  analysis  as  for  the  temperature  equilibration  term^  discussed 

in  Sec.  2,  because  there  is  no  simple  analytic  solution  to  Eq.  (23). 

A straightforward  but  sometimes  very  slow  running  approach  would  consist 

in  splitting  the  hydro  time  step  for  that  calculation,  recalculating 

P , at  each  of  this  fractional  time  step  and  reaching  small  temperatures 
rad 

in  a smooth  fashion,  as  shown  in  Figure  3h.  However,  a much  faster  and 

still  adequate  way  is  to  truncate  the  energy  loss  due  to  when 

reaches  a small  value,  at  time  interval  ‘t*.  We  then  set  T to  that 

e 

value  and  take  P . ^t*  as  the  energy  radiated  in  the  time  interval  ?t. 
rad 

This  approach  is  shown  in  Figure  3c.  Besides  being  very  fast,  it  has 
the  advantages  of  avoiding  negative  temperatures  and  also  of  not  over- 
estimating (in  fact,  underestimating  somewhat)  the  radiated  energy. 

Another  problem  where  the  atomic  physics  impacts  the  hydrodynamics 

is  the  relative  importance  of  the  ionization  energy  E , in  the  internal 

chem 

energy.  In  the  atomic  physics  part  of  the  code,  the  ionization  energy 
can  be  calculated.  From  that  number,  we  define  an  effective  adiabatic 
coefficient  v for  the  electron  thermal  plus  ionization  energy  by 


n kT 

e.  e e 


n kT 
e e 

(v-1) 


In  this  form,  the  electron  Internal  energy  represents  the  sum  of  the 
ionization  and  the  electron  thermal  energy.  When  energy  is  added  (or 
subtracted)  to  the  plasma,  a certain  fraction  is  taken  up  by  the  ioni- 
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zation  energy.  As  the  temperature  goes  up,  this  fraction  becomes 
smaller  and  smaller  (V  then  tends  toward  5/3).  At  low  temperatures, 

Echem  much  larger  than  the  thermal  energy  and  leads  to  a value 

of  V very  close  to  1.  This  in  turn  leads  to  problems  when  division  by 
(v-1)  occurs.  A minimum  v is  set  up  depending  on  the  element  in  order 
to  avoid  this  problem.  (Incidentally,  a correct  approach  to  this  problem 
would  be  to  use  a correct  equation  of  state  for  the  transition  state 
between  solid  and  plasma).  The  minimum  value  of  v is  determined  from  the 
ionization  energies  of  the  first  txro  levels.  For  example,  for  A;,  the 
minimum  value  of  Y used  in  the  code  is  1.15^. 

V.  CONCLUSION 

A new  wire  code  l-JHYRAC  has  been  developed.  This  code  provides  a 
nex^  tool  for  the  study  of  multiple-wire  arrays,  single  wires,  Z-pinches 
and  other  cylindrically  S3mimetric  plasma  configurations.  The  code  solves 
the  circuit  equations  with  the  plasma  included  in  a self-consistent  way. 
No  lumped  electrical  characteristics  are  used  for  the  plasma  in  those 
equations.  The  atomic/radiation  package  has  been  designed  to  Interact 
with  the  rest  of  the  code  in  a standardized  fashion  and  can  be  upgraded 
at  will,  depending  on  the  amount  of  computation  time  which  is  to  be 
spent  in  that  stage.  At  this  point  of  development  of  the  code,  a dis- 
tinction is  made  betx^een  a fast  "energetics"  package  and  an  elaborate 
diagnostic  package  x^hlch  can  be  used  as  a post-processor  for  data  gener- 
ated by  the  code.  As  a last  remark  (in  this  numerically  oriented  report), 
the  numerical  "fixes"  in  the  code  have  been  kept  to  a minimum  (and  exten- 
sively discussed  herein),  and  extensive  tests  have  been  made  in  order  to 
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choose  the  most  satisfactory  algorithm  to  deal  with  each  numerical 
problem  and  to  minimize  the  effect  of  these  choices  on  the  physics. 
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METAL  END  PLATE 

Fig.  2 — A schematic  of  the  wire  plasma  in  its  metal  return  cur- 
rent path.  The  current  loop  around  which  Faraday’s  law  is  inte- 
grated is  shown  in  dashed  lines. 
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-53) 


4 


EGSG,  INC. 

ALBUQUERQUE  DIVISION 
P.O.  BOX  10218 
ALBUQUERQUE,  NM  87114 

1 CY  ATTN:  TECHNICAL  LIBRARY 

FORD  AEROSPACE  S COMMUNICATIONS  CORP 
3939  FABIAN  WAY 
PALO  ALTO,  CA  94303 

(FORMERLY  AERONUTRONIC  FORD  CORPORATION) 
1 CY  ATTN:  DONALD  R.  MCMORROW  MS  G30 

1 CY  ATTN:  LIBRARY 

FORD  AEROSPACE  £ COMMUNICATIONS  OPERATIONS 
FORD  £ JAMBOREE  ROADS 
NEWPORT  BEACH,  CA  92663 

(FORMERLY  AERONUTRONIC  FORD  CORPORATION) 
1 CY  ATTN:  TECH  INFO  SECTION 

GENERAL  ELECTRIC  COMPANY 
SPACE  DIVISION 
VALLEY  FORGE  SPACE  CENTER 
GODDARD  BLVD  KING  OF  PRUSSIA 
P.O.  BOX  8555 
PHILADELPHIA,  PA  19101 

1 CY  ATTN:  JOSEPH  C-  PEDEN  VFSC,  RM  4230M 

GENERAL  ELECTRIC  COMPANY 
TEMPO-CENTER  FOR  ADVANCED  STUDIES 
816  STATE  STREET,  (P.O.  DRAWER  QQ) 

SANTA  BARBARA,  CA  93102 

1 CY  ATTN:  DASIAC 

INSTITUTE  FOR  DEFENSE  ANALYSES 
400  ARMY-NAVY  DRIVE 
ARLINGTON,  VA  22202 

1 CY  ATTN:  IDA  LIBRARIAN  RUTH  S.  SMITH 

ION  PHYSICS  CORPORATION 
SOUTH  BEFORD  STREET 
BURLINGTON,  MA  01803 

1 CY  ATTN:  H.  MILDE 

IRT  CORPORATION 
P.O.  BOX  81087 
SAN  DIEGO,  CA  92138 

1 CY  ATTN:  R.  L.  MERTZ 
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JAYCOR 

1401  CAMINO  DEL  MAR 
DEL  MAR,  CA  92014 

1 CY  ATTN:  ERIC  P.  WENAAS 

JAYCOR 

205  S.  WHITING  STREET,  SUITE  500 
ALEXANDRIA,  VA  22504 

1 CY  ATTN:  ROBERT  SULLIVAN 

KAMAN  SCIENCES  CORPORATION 
P.O.  BOX  7463 

COLORADO  SPRINGS,  CO  80933 


1 

CY 

ATTN: 

ALBERT 

P.  BRIDGES 

1 

CY 

ATTN: 

JOHN  R 

. HOFFMAN 

1 

CY 

ATTN: 

DONALD 

H.  BRYCE 

1 

CY 

ATTN: 

WALTER 

E.  WARE 

LOCKHEED  MISSILES  AND  SPACE  CO  INC 
3251  HANOVER  STREET 
PALO  ALTO,  CA  94504 

1 CY  ATTN:  LLOYD  F.  CHASE 

MAXWELL  LABORATORIES,  INC. 

9244  BALBOA  AVENUE 
SAN  DIEGO,  CA  92123 

1 CY  ATTN:  A.  RICHARD  MILLER 

1 CY  ATTN:  PETER  KORN 

1 CY  ATTN:  Wayne  Clark 

MCDONNELL  DOUGLAS  CORPORATION 
5301  BOLSA  AVENUE 
HUNTINGTON  BEACH,  CA  92647 

1 CY  ATTN:  STANLEY  SCHNEIDER 

MISSION  RESEARCH  CORPORATION 

755  STATE  STREET 

SANTA  BARBARA,  CA  93101 

1 CY  ATTN:  WILLIAM  C.  HART 

1 CY  ATTN:  CONRAD  L.  LONGMIRE 


MISSION  RESEARCH  CORPORAT I ON-SAN  DIEGO 
P.O.  BOX  1209 
LA  JOLLA,  CA  92038 

(VICTOR  A.  J.  VAN  LINT) 

1 CY  ATTN:  V.  A.  J.  VAN  LINT 
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NORTHROP  CORPORATION 

NORTHROP  RESEARCH  AND  TECHNOLOGY  CTR 
3401  WEST  BROADWAY 
HAWTHORNE,  CA  90250 

(DESIRES  ONLY  1 COPY  OF  CNWDI  MAT) 
1 CY  ATTN:  LIBRARY 

NORTHROP  CORPORATION 
ELECTRONIC  DIVISION 
2301  WEST  120TH  STREET 
HAWTHORNE,  CA  90250 

1 CY  ATTN:  VINCENT  R.  DEMARTINO 


PHYSICS  INTERNATIONAL  COMPANY 
2700  MERCED  STREET 


LEANDRO, 

CA 

94577 

1 

CY 

ATTN 

: DOC 

CON 

FOR 

BERNARD  H.  BERNSTEIN 

1 

CY 

ATTN 

: DOC 

CON 

FOR 

CHARLES  H.  STALLINGS 

1 

CY 

ATTN 

: DOC 

CON 

FOR 

PHILIP  W.  SPENCE 

1 

CY 

ATTN 

: DOC 

CON 

FOR 

Bemie  Lippmann 

1 

CY 

ATTN 

: DOC 

CON 

FOR 

SIDNEY  D-  PUTNAM 

PULSAR  ASSOCIATES,  INC. 

11491  SORRENTO  VALLEY  BLVO 
SAN  DIEGO,  CA  92121 

1 CY  ATTN:  CARLETON  H.  JONES  JR. 


R £ D ASSOCIATES 
P.O.  BOX  9695 
MARINA  DEL  REY,  CA 

1 CY  ATTN: 

1 CY  ATTN: 

1 CY  ATTN: 


90291 

C.  MACDONALD 
WILLIAM  R.  GRAHAM  JR. 
LEONARD  SCHLESSINGER 


SCIENCE  APPLICATIONS,  INC- 
P.O.  BOX  2351 
LA  JOLLA,  CA  92038 

1 CY  ATTN:  J.  ROBERT  BEYSTER 

SPIRE  CORPORATION 
P.O.  BOX  D 
BEDFORD,  MA  01730 

1 CY  ATTN:  ROGER  G.  LITTLE 


SRI  INTERNATIONAL 
333  RAVENSWOOD  AVENUE 
MENLO  PARK,  CA  94025 

1 CY  ATTN:  SETSUO  DDAIRIKI 
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SYSTEMS,  SCIENCE  AND  SOFTWARE,  INC. 

P.O.  BOX  4803 
HAYWARD,  CA  94540 

1 CY  ATTN:  DAVID  A.  MESKAN 

SYSTEMS,  SCIENCE  AND  SOFTWARE,  INC. 

P.O.  BOX  1620 
LA  JOLLA,  CA  92038 

1 CY  ATTN:  ANDREW  R.  WILSON 

TEXAS  TECH  UNIVERSITY 

P.O.  BOX  5404  NORTH  COLLEGE  STATION 

LUBBOCK,  TX  79417 

1 CY  ATTN:  TRAVIS  L-  SIMPSON 

TRW  DEFENSE  £ SPACE  SYS  GROUP 

ONE  SPACE  PARK 

REDONDO  BEACH,  CA  90278 

1 CY  ATTN:  TECH  INFO  CENT ER /S- 1 9 30 

VOUGHT  CORPORATION 
MICHIGAN  DIVISION 
38111  VAN  DYKE  ROAD 
STERLING  HEIGHTS,  MI  48077 

(FORMERLY  LTV  AEROSPACE  CORPORATION) 


1 CY 

ATTN: 

TECH  LIB 

NRL 

CODE 

2628 

- 20 

CYS 

NRL 

CODE 

6700 

- 1 

CY 

NRL 

CODE 

6750 

- 20 

CYS  (1  CY  CLASSIFIED) 

8 


